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Abstract 

In this paper we study the emergence of coherence in collective motion described 
by a system of interacting motiles endowed with an inner, adaptative, steering 
mechanism. By means of a nonlinear parametric coupling, the system elements 
are able to swing along the route to chaos. Thereby, each motile can display 
different types of behavior, i.e. from ordered to fully erratic motion, accordingly 
with its surrounding conditions. The appearance of patterns of collective motion 
is shown to be related to the emergence of interparticle synchronization and the 
degree of coherence of motion is quantified by means of a graph representation. 
The effects related to the density of particles and to interparticle distances are 
explored. It is shown that the higher degrees of coherence and group cohesion 
are attained when the system elements display a combination of ordered and 
chaotic behaviors, which emerges from a collective self-organization process. 
Keywords: Collective motion, swarming, coupled maps, synchronization, 
self-organization, chaos. 

1. Introduction 

Complex motion modes of collectives as a result of their constituent interact- 
ing entities occurs almost ubiquitously in nature and over the last decades it has 
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provided a common ground for cross-disciplinary investigations among Physics, 
Biology and Mathematics. The range of applications of such studies is indeed 
extensive . As a matter of fact, coherent patterns of collective 

motion found in distinct families of biological species such as fish schools, flocks 
of birds, swarms of insects and even colonies of bacteria 

have also been detected in granular matter systems, self-propelled particles with 
inelastic collisions and active Brownian particles in autonomous-motor groups 
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l2Ci |2l|, |22j. Early in the study of such collectives, 



modeling and simulation have been recognized as playing a crucial role in gain- 
ing insight of the mechanisms underlying such an emergence of global features 



from a set of simple rules [2J 



As the interest of the scientific community in addressing such kind of systems 
increases, minimal microscopic models have been recently introduced. Most of 
them consider systems of many interacting elements whose couplings, being 
in general nonlinear, can be either local or global. From a purely determinis- 
tic perspective, such systems can be represented as high dimensional dynam- 
ical systems, which can be either discrete or continuous in their time evolu- 
tion. This approach has been mainly spearheaded by Smale and collaborators 



HQ, 



211 ] . On the other hand, providing understanding of the connec- 
tion between macroscopic collective features and microscopic scale interactions 
in multi-particle systems is well within the principal aims of statistical physics. 
Therefore, naturally, models of self-propelled interacting particles have provided 
a fertile ground for study using concepts and tools of statistical physics. This 
kind of approach has been mainly adopted in models where randomness is intro- 
duced in the dynamics of the particles by means of a Langevin-type description 

Finally, another line of investigation aims at 
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providing descriptions and modeling in purely probabilistic terms where interac- 
tions obey probabilistic 'rules of engagement'. Notably, effective Fokker-Plank 
equations have been recently proposed for coarse-grained observables of 'agent 



systems' (see for example 12|, |31| and references therein for a more detailed 
presentation). 
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One of the earliest theoretical stochastic models of self-propelled interac ting 
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particles was introduced by Vicsek and collaborators as early as in 1995 
5j, which still possesses seminal value because of its minimal character. In 
Vicsek's model, point particles move at discrete time steps with fixed speed. At 
every time step, the different particles velocities are determined by the average 
of neighboring particles. In other words, Vicsek's model is an XY model in 
which the 'spins are actively moving' (see 32|). Furthermore, similarly as in 
ferromagnetic spin systems, Vicsek's model exhibits a phase transition as a 
function of both the particle density and the intensity of noise. For a detailed 
inve stig ation on the nature of such phase transitions we refer the reader to 



291 1151 130J. Further variations of Vicsek's model have been recently proposed 
to account for changing symmetries, adding cohesion or taking into account a 
surrounding fluid interacting with the particles 32j. 



Although the application of concepts stemming from statistical physics re- 
search has led to identify some universal properties existing in these classes of 



systems, such as s 
chronizing modes 
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symmetry breaking, phase transitions and syn- 



28l Ijj], the role of the individuals' internal dy- 



namics still remains veiled. 

Whilst the explicit consideration of inner control processes could increase the 
complexity of models of interacting motiles, it is a necessary conceptual step 
in developing further insights into the mechanisms underlying the emergence 
of coherence of motion in biological systems. Historically, models of group 
motion where particles adapt to their environment by means of an inner steering 
mechanism have been developed in the context of traffic modelling. For an 

n 

overiew of such models we refer to [33|. In the context of biological systems, 
inner states have been considered in order to model the emergence of coherent 
behavior in groups of fireflies 
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35j . More closely related to the problem here 



addressed is the study of the response and adaptation of populations of motiles 
to the information carried by external 'fields'. Recently, by assuming inner state 
dynamics, attempts in this direction have been reported in the study of bacterial 
chcmotaxis 36 



37] and biologically inspired collective robotics 
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In the present work we introduce a purely deterministic model where, in 
analogy to the Vicsek's class of models, particles can display phenomenological 
random-like motion and exhibit a sharp change of coherence of motion as a 
function of the particle density. Furthermore, in contradistinction to Vicsek's 
class of models, no boundary conditions are considered, since a feature of the 
group's cohesion is that it is built by the collective dynamics 'per se'. Even in 
the absence of explicit interparticle attraction, it is the coordinated effectively 
synchronized collective motion that keeps the group together. In our model ev- 
ery motile is endowed with an inner 'steering' variable that evolves according to 
an heuristic discrete-time equation. For the sake of simplicity, such an evolution 
law has the structure of the logistic map. The latter provides a suitable, well 
understood, combination of chaotic and ordered behavior to account for the 
coherence and novelty aspects observed in real collectives of motiles. Commu- 
nication between a particle and its environment occurs via a control parameter 
that tunes its value according to the external states of the surrounding particles. 
At the microscopic level the features of motiles are summarized by the following 
conditions, along the lines of [3]: 

(a.) Each element has a time dependent internal state and spatial position. 
((3) Each element is 'active' in the sense that its internal state can exhibit chaotic 
behavior, both in presence and absence of interactions with other particles. 

(7) The dynamics of the internal state of a given element is determined by local, 
short range interactions, effectuated within a neighborhood of a characteristic 
radius. 

(8) The interparticle interaction depends on the internal states of the partici- 
pating particles. 

These general rules have been found to give rise to nontrivial emergent behavior 
which can not be readily deduced from the microscopic parameters of the system. 

At the collective or 'macroscopic' level, the characteristic emergent phenom- 
ena observed in such kind of locally coupled systems are mainly described by 
the notion of 'clustering', either in real or in state space. As it has been reported 
in [nl , distinct classes of clustering behavior accompany, in a generic way, such 



4 



a coupling: 

(i) Elements forming a cluster merge in and out of the cluster. 

(ii) Elements can remain separated from neighboring clusters but they form a 
bridge between distinct clusters facilitating information flow, exchange of ele- 
ments between clusters and adding cohesion. 

(iii) Presence of independent clusters separated by distances larger than the 
interaction, with elements rarely merging in and out of the clusters amidst 
them. 

(iv) Cluster - cluster interactions such as aggregation, segregation and compet- 
itive growth between various sized clusters. 

In this work, we shall focus on the description and quantification of emergent 
collective motion based on a clustering index that accounts simultaneously for 
both, the degree of spatial clustering and the degree of interparticle alignment 
of velocities. 

The paper is organized as follows: In Section [2J we introduce a general for- 
mulation of the model. Next, in Section [3] we address the possible types of 
stationary behavior in the motion of an individual particle, as well as the basic 
phase-locking synchronization process that results from local interparticle pair- 
interactions. Section[4]presents the case of the many particle system. Its typical 
evolution patterns, regimes of motion, degree of synchronization and the depen- 
dency on both, density of particles and interparticle distance, are investigated. 
Finally, Section [5] concludes the present work with a brief summary, discussion 
of the results and possible further extentions. 

2. Formulation of the model 

We consider N particles, labeled through an index i = 1,2, . . . ,N, whose 
positions at time t are denoted by the vectors {rj}. They evolve on a plane (two 
dimensional motion) where their positions change simultaneously at discrete 
time steps At, according to 

f* +A4 = + At. (1) 
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Similarly as in Vicsek's original model [29l.l5|. we assume here that at every time 
step the speeds of all particles are equal to a common constant value 



Changes in the particles' velocities occur via an inherent steering mechanism 
which can be expressed in terms of a two dimensional rotation matrix T\ as 



Assuming that the motion of each particle is governed by an inner steering 
process, we endow each particle i with a variable 9\ determining, at every time 
step, the phase of the rotation matrix T\= T(0£). Phases 6\ are assumed to 
take values in an interval [— Ao, Ao], where the maximum rotation angle Ao is 
taken to be a small fraction of n. Furthermore, let us consider the evolution of 
the rotation phase, for each of the particles, as determined by an equation of 
the general form 



Here, the function &\ is introduced to model the response of the particle i to 
the influence exerted by its 'environment', i.e. the set of positions and velocities 
of all the other surrounding particles. In the present framework we require the 
functions <&\ to fulfill the following generally admitted conditions: 

a. Two particles will interact provided they are inside a neighborhood of fixed 

radius R. 

b. The intensity of the interaction between two neighboring particles should 

decay as a function of their interdistance. 

c. In absence of any neighbor, particles should follow an unbiased, completely 

erratic trajectory. 

d. Frontal collisions between pairs of particles should be hindered. 

e. The interactions within a group of particles should lead to the emergence of 

coherent patterns of collective motion. 




(2) 




(3) 
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/. The cluster formations made by particles should maintain a certain degree 
of cohesion. 
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Figure 1: a) Particles interact via short range, local interactions extending over a neighborhood 
of radius R. b) The graph of the weight function {8| where d = \\r\ — r\ || denotes the 
interparticle distance. 

The coalescence of coherence and novelty observed in real collectives of 
motiles suggests to consider an heuristic function <fr\ being able to display be- 
haviors ranging from order to fully developed chaos. As well known, the logistic 
map Xt+\ = fJiXt(l — Xt) with x S [0,1] and < /i < 4 [3], is a minimal 
representative model for the class of discrete dynamical systems unfolding the 
route to chaos via a period-doubling bifurcation cascade. Hence, we propose 
the following function 




#i = Ao-2.Ao-#.^l-^J ) (.-,, 

which is obtained by introducing the change of variable 

xt -> \ U - and /i -> Acj> (6) 

in the logistic map. Here, the functions 4>\{^ |rf ,vj j ^ in ([5J should embody 
the coupling between particles in such a way that conditions a - f are satisfied. 
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As it is shown throughout the forthcoming sections, the latter is achieved by 
the following function 

4=1 , (7) 

I 1 if n» = 

where n$ is the number of neighboring particles counted within a neighborhood 
Di{R) of fixed radius R around each particle i. Such neighborhood guarantees 
the validity of condition a. The interaction of neighboring particles i and j is 
weighted by a suitably chosen function Wi_j which aims at satisfying condition 
b. In particular, we shall assume that Wij is given by 

R K 

w i,j = : ~k- (8) 

R K + (\i?t-r*\\) 

Here, the parameter K controls the degree of dependence of the coupling on the 
interparticle distance. Notice that function (|5J) approaches a step function as K 
increases (see Fig. H{a)). As it turns out, for K ^> 1 the effect of different inter- 
particle distances practically disappears, as the contributions of all interparticle 
couplings, within the neighborhood, tend to be equally weighted. 

3. Local Dynamics 

We proceed now to the more 'microscopic' local level in order to elucidate 
the underlying dynamics of particle-particle interactions. It is instructive to 
consider one particle in isolation and subsequently a single pair of particles 
and their interaction. In a sense this is analogous to the statistical mechanical 
treatment of 'dilute gas' where the particles' collisions, being rare, are described 
very accurately by their binary collisions. Since, at very low densities, binary 
interactions are the most dominant contributions in the present model, we shall 
consider such a scenario in this section. 

3.1. The one particle case 

Let us consider first a single particle with its parameter <f>\ = <fi being a 
constant. In the present case, the whole information about the possible types 
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of particle trajectories is contained in the bifurcation scenario for the variable 
9t as a function of <f>. Such bifurcation diagram is qualitatively the same as 
in the logistic map, as it is carried through by the transformations © and 
<j> = /z/4. As shown in Fig. [5^), the bifurcation diagram of the rotation phase 
depicts fixed points in the interval [0, |], the emergence of multiple period orbits 
[MPO] in the interval (|, F p ] (where F p « 0.8925 corresponds to the well-known 
Feigenbaum accumulation point |3|]), onset of chaos at F p , coexistence of weak 
chaos and periodic motion in (F p , 1) and finally, fully developed chaos at <f> = 1. 
It is emphasized here that for < <p < \ t he rotation phase attains a maximum 
value corresponding to the trivial stationary solution 9 = A of eqs. (@| - ([5]). 
For j < <p < | the rotation phase is given by 9 = Ao- The particle 
thereby follows closed trajectories with a radius of curvature growing as <\> — > i. 
At <j> = I a change in the sign of the rotation occurs, from anticlockwise for 
< 4> < h t° clockwise for | < <j> < |. In order to characterize the trajectories 
of the particles in the interval [|, 1], it is worth considering the mean rotation 
phase 

W = -£* (fc) 0) 

/c=i 

where m denotes the period of a certain orbit and 9^ denotes its k^ 1 element. 
In the (^-interval (|, 0.862], which corresponds to orbits of period 2, the mean 
rotation phase reads (9) = Providing an analytic expression for (6) in the 
case of orbits of period higher than 2 is a cumbersome task and therefore, for 
0.862 < (f> < 1, we calculate numerically the values of (9). These results are 
presented in panel (b) of Fig. [2] For <fi values in the interval (0.862, F p ), the 
trajectories are closed and clockwise. For <f> m the interval (F p ,l), the mean 
phase (9) displays highly irregular behavior, where the alternation of windows 
of weak chaos and periodic behavior gives rise to a fine structure profile. In 
the case of weak chaos, the particle performs a biased chaotic walk leading to 
quasi-circular trajectories. All along the interval (Fp, 1), the particle motion is 
clockwise on the average. The peaks observed in Fig. Hfb), where (9) attains 
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Figure 2: a) Bifurcation diagram of 0t/Ao as a function of the parameter (f>. b) Similar plot for 
the mean phase (<9t)/Ao (see eq. 10). The inset of panel (b) provides a zoom-in for {f?t)/Ao 
in the ^-interval (0.89, 1]. 

higher values, are indicative of the presence of windows of periodicity. Finally, 
at the limiting case of fully developed chaos (i.e. <j> — 1), the rotation phase 
takes values distributed within the interval [— A , A ], according to the invariant 
density 

g,(M = i) = * fl2 

TTy/A 2 - e 2 

of map (|4]) - (|5]) . The even character of g s entails that the particle describes a 
strongly chaotic walk with (8) =0, as it is confirmed by Fig. [2jb). This way 
condition c above is satisfied since, according to (0, all isolated individuals 
evolve with <f>\ = 1. This analysis has been carried out for a single particle with 
an arbitrary fixed value of (f>. Yet, the same categorization readily applies to 
the case of many interacting particles, as long as their coupling parameters 
converge towards a quasi-stationary value. The main results obtained from the 
one particle case are summarized in Table I. 
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Table 1: Classification of the asymptotic behavior of 9 and of the 
particle trajectories, for characteristic domains of the parameter (f>. 



INTERVAL: PARAMETER RANGE 


ASYMPTOTIC BEHAVIOR OF 6 


TYPE OF 
TRAJECTORY 






closed, 






anticlockwise 






trajectories for 






< 4> < 1/2; 


A: 6 [0,|) 


stationary state 


straight 
trajectories at 
<\> — 1/2; clockwise, 
closed trajectories 
for 1/2 < (f>< 3/4 


B:0 6[|,F P ) 


multiple period orbits 


closed, clockwise 
trajectories 






biased chaotic 






trajectories 


C: 6 [F p , 1) 


weak chaos 


(quasi-closed, 

clockwise 
trajectories) 


D: (f> = 1 


fully developed chaos 


unbiased chaotic 
trajectories 



3.2. The two particle case: interparticle synchronization 

The case of two particles is useful to consider in order to illustrate some 
of the control features exercised by the coupling function ([7]) which, in this 
particular case, is the same for both particles, <jyp — <f>t = (f^ 1 ' 2 ^ ■ Unless 
stated otherwise, the set of parameters we use throughout the rest of the paper 
are: s = 1, At = 1, R = 500, A = ±tt and K = 1. 

First, it is straightforward to show how condition d is satisfied. Figure [3] 
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depicts the plot of the trajectories of two particles initially being out of the 
interaction zone and approaching each other along the same axis. Under such 
conditions, when they both enter the interaction zone, the difference of velocities 
approaches its maximum value \\Vf^ — II — 2s and Wij — i. Therefore, as 

(1 2) 

soon as both particles start interacting, the parameter (f> t ' drops sharply from 
unity to |, as it is shown in Fig. EJb). According to Fig. Ufa), a value of 
< I entails a single fixed point for both phases 6^ and 9^ . As a result, 
the difference l^ 1 -* — will tend to vanish after a brief transient. The latter 
can be realized by comparing Figs. [3jb,c,d). 
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Figure 3: a) Trajectories of two particles initially approaching each other frontally. b) Plot of 
</»j 1,2 ' versus time t, corresponding to particles of panel a), c) Similar plot for 6^ /Aq and 
0^ 2 ' /Ao. The parameters used are s = 1, At = 1, R = 500, Ao = ^jtt and K = 1. 

Moreover, since for Hrf'-rf^l -» the weight W(i.2) — >■ 1, a further decrease 

(1 2) 

in cj)\ occurs as both particles get closer (see the decrease towards a minimum 
value in the plot of Fig. Efb)). As it was discussed in Subsection GHU a low value 
of <^ < \ leads both phases to reach their maximum value Aq (see Figs. [3^) 
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and b)). This fact is demonstrated by the onset of the plateau appearing in 
Figs. [2Ic,d). Once the phases of the particles reach their maximum value, 
simultaneously with their minimum phase difference, they follow their course 
away from each other. As it turns out, condition d is satisfied as both particles 
turn away sharply instead of colliding frontally (see Fig. [2a)). Eventually, both 
particles leave the interaction zone and display again fully erratic trajectories. 
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Time 



Figure 4: a) Trajectories of two particles located initially inside the interaction zone and 
approaching each other with a moderate difference of initial velocities. In panel b) we present 
the plot of (pi ' ' versus time t corresponding to particles of panel a). In panels c) and d), 
we show similar plots for 6^ j Ao and 0^ / Aq respectively. The parameters we used in these 
plots are the same as in Fig. [3] 

Let us consider now the case of two particles approaching each other inside 
the interaction zone with a small difference in their initial directions. Figure 
[4ja) shows the trajectories of two particles in such a typical situation, while Fig. 
0Jb) depicts the evolution of the coupling parameter 4>\ ' ' and Figs. |U[c,d) show 
the evolution of the rotation phases 9^ and 9^ of both particles. Initially, a 
sharp decrease in the difference of directions \\v^^ — vf^ 1 1| occurs, which appears 
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in the plot of Fig. QJb) as an initial abrupt increase of <f>\ ' . Since (jr t ' attains 
a value higher than | , the rotation phases of both particles subsequently enter 
into the MPO and chaos regimes. The ongoing evolution process is driven by 
two counteracting tendencies: The first one is a trend towards strong synchro- 

(1 2) 

nization if <j)\ decreases below Fp. The second one is a trend to spread away 
if an alignment occurs, i.e. if (p\ — > 1. The presence of these two trends is 
revealed when comparing panels (a), (b) and (c) of Fig. @] the values of 0j and 

(2) • (12) 

9 t tend to further localize when (jr t : decreases, while bursts of chaos appear 

f 1 2) 

whenever <p). ' ' « 1 (corresponding to an alignment). As it is depicted in Figs. 
0Ja-c), after a transient self-organization process, the evolution of both particles 
undergoes a transition towards a highly coherent dynamic regime. The initially 

(1 2) 

irregular behavior observed in the course of the evolution of <fi t suddenly dis- 
appears and, instead, oscillations emerge within the ^-interval (0.85,0.9). Such 
oscillations of ^j 1 ' 2 ^ induce on both phases 8^ and 8^ a regular alternation 
between windows of high and low period orbits (see Figs. [4] (c,d)). Once such 
a new regime is attained, both particles tend to move in perfect synchrony and 
their phase difference \9^ — 9^\ equals zero, as it is readily shown in Fig. [5j 
As it turns out, for the case of two particles, it is clear that at longer term 
the interplay between both tendencies (ordering and chaotic phase spreading) 
constitute a 'control mechanism' that underlies the fulfillment of conditions e 
and /. Using the set of parameters as in Fig. |3]for different settings of initial 
conditions, we found that phase synchronization, as it is observed in Fig. [5l 
occurs provided the initial conditions are such that both particles never leave 
the interaction zone. In Section|4l we show that phase synchronization is robust 
in the case of many particles as well. 

4. Emergent global dynamics and patterns of motion 

Having seen the role of spontaneous synchronization in binary interactions 
we can now consider the full collective motion problem. As we shall demonstrate 
further in this section, synchronization is an underlying mechanism for coherent 
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Figure 5: Corresponding to phases 9^ and 6^ of Fig. 2] the plot of the time evolution of 
the phase difference (0 t (1) - t (2) )/A o . 

and cohesive motion. Another feature of this intrinsic and coherent steering, 
which keeps a large fraction of the group together, is that the role of boundary 
conditions becomes really secondary. Furthermore, we emphasize the fact that 
the aggregates of motiles keep its coherent motion even in absence of boundaries. 

So, let us address now the case of N particles allowed to move on the infinite 
plane. We consider, thereof, that the initial positions and velocities of the par- 
ticles are uniformly randomly distributed within a square of side Rq. Similarly, 
the rotation phases are distributed in the interval [— Ao, Aq]. 

In order to capture the spatio-temporal features of the system, such as per- 
sistent and/or flickering clustering as well as patterns of cohesion around tempo- 
rary foci, it is of interest to introduce an index that takes into account both the 
spatial relations between neighboring particles and the degree of velocity align- 
ments. To this aim, it is instrumental to visualize the system as an undirected 
graph, where every particle is its node and the edges of the graph are established 
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according to their intcrparticle distance and to the degree of alignment between 
particle velocities. 

In particular, we consider a pair of nodes as connected if, in addition of 
being neighbors, the inner angle between their velocities is less than Ao. Under 
these assumptions, the time dependent adjacency matrix M t associated to such 
a graph reads 




^ - II < R and \\v\ - 4 1| < sA 



(10) 



In terms of this matrix, one can quantify the degree of collective alignment 
resulting from interparticle interactions by means of the following index 

(11) 

J- M 

where T M = ±N(N - 1) is the maximum number of possible connections and 
T t = \Tr [M t 2 ] (with Tr denoting the matrix trace and M t 2 the square of the 
adjacency matrix) gives the number of actual connections at time t. Hereafter 
we refer to index a t as the Alignment Clustering Index (ACI) at time t. Clearly, 
ACI is an indicator of the degree of coherence of motion of the group of particles. 
Furthermore, it provides information on the degree of clustering of particles 
in space and thus it constitutes an indicator of the degree of group cohesion. 
For instance, strong oscillations in the evolution of the ACI indicate a poorer 
group cohesion, since its fluctuations reveal that a great amount of particles are 
merging in as well as escaping from cluster formations. 

4-1- Typical patterns and regimes of motion 

In order to demonstrate some typical patterns of characteristic behavior 

of our model, let us consider a collective of particles initially distributed in a 

squared region of side R = —R. Unless explicitly mentioned, we shall consider 

3 

the number of particles N = 300 and the rest of the parameters to be the same 
as in Fig. [3] 

Figure [6] shows the plot of the positions and velocities of the set of particles 
at different times in a single realization. Corresponding to the system of Fig. [51 
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Figure 6: For a system of TV = 300 particles, in absence of boundaries and initially uniformly 
randomly settled within a squared region of side Ro = , plots of the positions and velocities 
at a) t = 0, b) t = 2 X 10 3 , c) t = 6 X 10 3 and d) t = 1.5 X 10 4 . The parameters are the same 
as in Fig. [3] 
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panel (a) of Fig. [7] shows the evolution of the ACI and panel (b) the evolution of 
the first and second moments of the coupling functions ([7]) of all particles. One 
can distinguish between three qualitatively different stages of collective motion, 
which can be characterized by means of the ACI behavior: 

1. Starting from a disordered state, Fig. Eta), the system undergoes a strong 
self-organization process, from t = 0uptot«2x 10 3 . This is char- 
acterized by low ACI values (see Fig. [JJJa)), and by large oscillations at 
the level of the coupling functions, Fig [7th). At the end of this stage the 
collective motion becomes further organized and the group of particles 
shapes a circular rotational pattern (see Fig. Etb)). 

2. The second stage, which extends approximately from t w 2 x 10 3 up to 
t « 1 x 10 4 , is mainly characterized by a marked increase of the ACI 
(see Fig. UJa)). The coupling functions tend to stabilize around a fixed 
value, as exhibited by the decrease in the amplitude of oscillations of the 
mean coupling (see Fig. [Tfb)) . Along this stage, an increasing number 
of groups of highly aligned particles are formed (see Fig. Etc)), thereby 
fulfilling condition e. 

3. At the final stage, the system attains a quasi-stationary regime, at t w 
1 x 10 4 , along which the ACI stabilizes and the coupling functions remain 
nearly constant (see Figs. Eta,b)). In such a regime, the standard de- 
viation of the average also tends to stabilize around a rather low value, 
thereby indicating an enhancement in the coherence of the behavior of the 
particles (see Figs. [T^c) andEJd)). Clearly, group cohesion is evidenced 
by the fact that high values of the ACI are sustained for times as long 
as t = 1 x 10 4 and beyond. As it turns out, condition / is also satisfied. 
It is a remarkable fact that the self-organization process leading to such 
state of higher coherence drives also the average of the coupling functions 
towards the region around the Feigenbaum point (which is represented by 
a dashed line in Fig. [JJJb)). In other words, a higher degree of coherence is 
achieved by a process that exhibits a combination of weak chaos and fur- 
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ther ordered behavior. This question will be discussed later in Subsection 
14.31 where we shall address the cohesion features of the presented system. 




Time 

Figure 7: Corresponding to the realization of Fig. [6] the evolution of the ACI is shown in 
panel a). In panel b), the mean value of the coupling functions (j>\ of all particles (where 
bars provide the corresponding standard deviation) is presented. The dashed line of panel a) 
corresponds to the value of the Feigenbaum point. 

4-2. Synchronization processes in the collective 

In Subsection 13.21 we have shown that the emergence of coherent behavior in 
the case of two particles arises from synchronization at the level of their rotation 
phases. Thus, it is natural to inquire on the degree of phase synchronization 
underlying the coherent aspect of the dynamics of the iV-particle system as it 
is exhibited in the plots of Figs. |5]and[7] In order to assess interparticle phase 
synchronization in a quantitative manner, we shall consider the ensemble of the 
distinct interparticle phase differences A0y = 9i — 9j with i ^ j. We shall, thus, 
monitor the evolution of the system in terms of its synchronized population 



19 



fraction. Such task is achieved by quantifying the percentage Pt(A9) of pairs 
of particles with a given phase difference AO at time t. Thereby, assessing 
the synchronization process amounts to follow the evolution of Pt(AO) on the 
A#-intcrval [— 2Ao,2Ao] starting from an initial distribution Pq. Clearly, full 
synchronization will be present whenever Pq converges asymptotically towards 
a delta distribution at A8 = 0. 




Phase difference 



Figure 8: Plot of the evolution of two different initial distributions of phase differences for 
iV = 300 particles. Upper panels: Having all particles the same initial phase 8, distribution 
of phase differences at (a) t = 0, (b) t = 2 X 10 3 , (c) t = 1 X 10 4 , (d) t = 4 X 10 4 and (e) 
t = 1 X 10 7 . Lower panels: for particles with initial phases uniformly randomly distributed, 
at (f) t = 0, (g) t = 2 X 10 3 , (h) t = 1 X 10 4 , (i) t = 4 X 10 4 and (j) t = 1 X 10 7 . The 
parameters used in this figure are the same as in Fig. \3\ Such numerical experiments supply 
ample evidence of the existence of an attracting distribution. 

Figure [8] exhibits the plot of the distribution P t (A9) at different times (from 
t = up to t = 1 x 10 7 ), starting from two different initial phase distributions: 
in upper panels (from (a) to (e)), the evolution of an initial delta distribution. 
In lower panels (from (f) to (j)), similar plots corresponding to an initial distri- 
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bution obtained by assigning uniformly distributed random values to the initial 
set of phases {6\} (the triangular shape in the initial distribution in Fig. [SJf) fol- 
lows from the correlations introduced by taking the differences between phases) . 
After a very short transient period, the initial correlations of phases correspond- 
ing to both distributions are destroyed and new ones are built as a consequence 
of the initial strong self-organization process described in Subsection 14.11 As it 
is shown in the upper and lower panels of Fig. [3J phase synchronization grad- 
ually emerges. At t = 1 x 10 4 the percentage of pairs with phase difference 
|A0| < Ao/10 is 38% for the case of the initial delta distribution (Fig. M,c)) 
and 40% for the case of the homogeneous initial distribution (Fig. [SJh)). As 
the process evolves, the interparticle phase synchronization is further enhanced. 
For instance, in the case of the delta initial distribution, the percentage of pairs 
with |A0| < Ao/100 at t = 1 x 10 7 equals 36% (Fig. [Sfe)) and it equals 60% 
for the homogeneous one (Fig. [SJj)). Since similar results are obtained for very 
different types of initial distributions, one concludes that the synchronization 
phenomenon observed in Fig. [8] is robust for parameter values as in Fig. [3j 
However, as it is shown in Subsection |4~3j similar self-organization processes are 
observed in a wide range of if-parameter values. 

Regarding the relation between the emergence of synchronization and the 
gradual formation of clusters of aligned particles in space (see Figs. [U[c,d)), it is 
interesting to inspect the evolution of the relation between pair synchronization 
and interparticle distance. Figure El depicts the plots of phase difference vs 
interparticle distances for the particle pairs accounted for in the histograms of 
panels (f) - (j) of Fig. [H Initially, all pairs are close to each other, having 
phase differences randomly distributed in the interval [— 2Ao,2Ao] (Fig. Eta)). 
At t = 2 x 10 3 , particles are further spread away and higher density 'clouds' 
appear within a narrower A0-interval (Fig. EKb)). The evolution towards a 
higher degree of organization leads to cluster formation in both space and phase- 
difference. As it is displayed in Fig. E£c), at t = 1 X 10 4 most of the pair particles 
are distributed in two large clusters. Strong synchronization (A0 « 0) occurs 
between particles within the same cluster, as well as between particles belonging 
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to different clusters. Also, many pairs inside and between clusters exhibit phase 
differences around specific values which correspond to the secondary peaks in 
the distribution of Fig. [5Jh). For a longer time, i.e. t — 4 x 10 4 , the main 
clusters split into smaller lumps within which particles are further synchronized 
with each other (Fig. EJe)). At t = 1 X 10 7 , a single compact cluster of fully 
synchronized pairs ( A9 = 0) still remains. As it turns out, the presence of strong 
synchronization at long term entails also a high degree of group cohesion. 

In comparison with other models of collective behavior, such as Vicsek's, one 
naturally expects the degree of group cohesion to decrease along with the initial 
density of particles. This question, as well as the dependence of the cohesion 
on the parameter K in the weight function (|H]), are addressed in the following 
Subsection. 




d 



Figure 9: Corresponding to plots f) to j) of Fig. [8] plot of phase differences (A(9/Ao) versus 
interparticle distances at a) t = 0, b) t = 2 X fO 3 , c) t = 1 X 10 4 , d) t = 4 X 10 4 and c) 
i = lx 10 7 . 
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Jj.,3. Group cohesion: the role of density and interparticle distances 

In the context of animal societies, sharp transitions of collective behavior 
have been reported to occur at specific values of the population density (see 



for instance 



18, 



y, Q, Q B, I 
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22]). As it has been already mentioned 



in the introduction, models belonging to the same class as Vicsek's exhibit an 
order-disordered phase transition at a critical value of the particle density. In 
the model introduced here, the radius of interaction R imposes a characteristic 
density p* = N/R 2 . Thus, it is natural to study the degree of coherence and 
cohesion in a group of particles for density values above and below p*. Figure 
[10] depicts the plot of the ensemble average of the ACI for a group of N = 100 
particles at t = 4 x 10 4 and for different values of the ratio p/p*. This plot 
exhibits a steep monotonic ascent of the averaged ACI occurring for values 
p/p* ~ O(10 -2 ), followed by the onset of a plateau extending for values p/p* > 
1. Such a behavior is best fitted by a logistic function of the form 

b 

Wp/p*)) t= 4 XW * = a + cp - (12) 

where a = 0.170, b Q = 0.003, c = 0.173 , b = (b - a)c p and p = 2.460 (as it is 
represented in Fig. [TU]by a continuous gray line) with a correlation coefficient 
of about 0.976. 

It must be noted here that the choice of the logarithmic scale for the ratio 
p/p* serves in elucidating the abrupt change observed. The change of behavior 
exhibited in Fig. [TQl as described by the logistic form of Eq. (fT2]l . differs in 



nature from the phase transition reported for Vicsek's class of models [29J, 132 1. 
In addition, we stress the fact that the abrupt change in the degree of coherence 
shown in Fig. [TU] arises naturally from the interparticle interactions, through 
concerted self-adaptation in the steering parameters of the particles. 

On the other hand, the cohesion of the group is controlled additionally by 
the parameter K in the weight function ((8|) . Figure [TTJ shows the plot of the 
ensemble average of the ACI at t = 4 x 10 4 for different values of K . One 
observes that stronger cohesion, as gauged by the ACI, is attained for small 
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p/p* 



Figure 10: Plot of the ensemble average of the ACI, (at(p/p*)}, of N = 100 particles at 
time t = 4 X 10 4 (dots) for different values of the logarithm of the ratio p/p*. Here the 
average was taken over an ensemble of 200 uniformly random initial conditions. The rest of 
the parameters are the same as in Fig. \3\ Bars in the plot indicate the standard deviation. 
The curve described by these points is best fitted by a logistic function in the form of Eq. 

G3. 
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Figure 11: For a group of N = 100 particles, plot of the ensemble average of the ACI, at 
t = 4 X 10 4 for different values of the parameter K. Similarly as in Fig. 1101 the average is 
computed over an ensemble of 200 uniformly randomly distributed initial conditions where 
error bars indicate one standard deviation. 

values of K. For < K < 1.25, the average ACI takes values around 15 — 20%. 
A decrease in the averaged ACI is observed in the interval < K < 5 which 
is followed by an increase for 5 < K < 10. For K values greater than 10, the 
averaged ACI gradually decreases. It is noticeable that the change in the overall 
behavior of the ACI around K = 5, from decrease to increase, coincides with 
the change of sign in the curvature of the graph of weight functions Wi j of 
([U when varying K. Since the form of the weight functions ultimately affects 
the behavior of the coupling functions <fii of ([7]), it is worth considering the 
ensemble average of the histograms of the coupling functions of all particles. In 
particular, we shall focus on the main regions described in Table I, namely A: a 
steady state, B: multi-periodicity, C: weak chaos and D: fully developed chaos. 
Figure [T^] displays the ensemble average of the histograms obtained at dif- 
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ferent times and for different vaiues of K. Initially (first column of Fig. [T2")) , 
regardless of the value of K, one observes that the coupling functions take values 
in the region where a steady state exists. For K — 1 (first row of Fig. [T2")) . the 
values of the couplings gradually evolve towards the region of multi-periodicity 
and beyond. At t — 4 x 10 4 one observes that the bulk of particles is equally dis- 
tributed between the multi-periodicity and weak chaos regions. However, 10% 
of the particles also appear in the region of fully developed chaos, which reveals 
particle escapes away from the center of mass. In the cases K = 5, K — 10 and 
K — 15 (rows 3, 4 of Fig. [P2)) , at long times, more than 20% of the particles have 
escaped and those that remain confined exhibit ^-values mostly in the regions 
corresponding to stationarity and multi-periodicity. Summarizing, in the case 
K = 1, a higher degree of coherence and group cohesion is attained through 
a process encompassing balance between weak chaos and ordered behavior. In 
the cases K = 5, K = 10 and K = 15, where most of the interacting particles 
display further ordered behavior, the degree of cohesion is poorer. In particular, 
in the case K = 5 almost 40% of the particles have escaped at t — 4 x 10 4 (as 
it is indicated by an increased number of particles with <pj 6 D). Consequently, 
for K — 5, the averaged ACI attains a very low value of approximately 0.05 (see 
Fig. HT|) . For 5 < K < 10, a further rigid steering motion arises, as an increased 
number of particles having <j)j 6 A tend to display short concentric trajectories. 
Thereby, the number of escapes is smaller and the average ACI becomes higher 
than in the case of K = 5 (see Fig. UTI) . The latter provides an insight into the 
origin of the minimum of the averaged ACI observed around K = 5. 

5. Conclusions 

In this work we have introduced a minimal model of motile particles which 
takes into account an intrinsic steering mechanism. Its main feature is that it 
exhibits the emergence of patterns of coherent collective motion. Since each 
particle is endowed with an intrinsic mechanism, it allows it to adjust its tra- 
jectory according to the surrounding conditions. With this extra feature the 
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Figure 12: Plots of the ensemble average of histograms obtained from a group of N = 100 
particles for the ^-intervals of Table I which correspond to A: a steady state, B: multi- 
periodicity, C: weak chaos and D: fully developed chaos (the percentage of particles in D 
indicates the proportion of isolated particles). The average was taken over an ensemble of 200 
initial conditions. Columns from left to right in this plot present results obtained for t = 1, 
t = 2x 10 3 , t = 1 X 10 4 and t = 2 X 10 4 while rows from top to bottom correspond to results 
obtained for K = 1, K = 5, K = 10 and K = 15. The rest of the parameters are the same as 
in Fig. O 

model can be of utility as an augmented inert particle model. The adaptation is 
achieved by changes which are determined by a map of the logistic-map family, 
which displays fully developed chaos in absence of interparticle interactions. As 
it turns out, isolated particles behave as chaotic walkers. On the other hand, 
particles within a neighborhood of fixed radius interact with each other by es- 
tablishing nonlinear couplings between their corresponding logistic maps. The 
coupling between a particle and its surrounding neighbors is embodied by a 
function that depends explicitly on the positions and velocities of the set of 
involved particles. The coupling functions play their role in self- adapting the 
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controlling parameters of the logistic map at the individual level of each motile 
by 'tuning in' the behavior of each particle. This behaviour can range from sin- 
gle to multiple periodic and chaotic regimes. The explicit form of the coupling 
functions is such that frontal collisions arc hindered and that neighboring par- 
ticles undergo a self-organization process leading to the emergence of coherent 
collective motion. The whole process is characterized by an index, called ACI 
(and denoted by a herein). This graph index, captures the essential features of 
the system's time evolution, since it quantifies the degree of cohesion and align- 
ment of particles. It is shown that the self-organization process of steering leads 
the collective of these motiles towards a regime of coherence and group cohesion. 
The latter is shown to be essentially related to the phase synchronization that 
emerges between the inner steering mechanisms of the particles. 

Also, the degree of group coherence is studied as a function of a parameter 
K. This parameter sets the dependence of the interactions on the interparticle 
distances. Such a study shows that higher levels of coherence and group cohesion 
occur in cases where most of the particles exhibit a balanced combination of 
ordered motion (multi-periodicity) and weakly chaotic behavior. 

Similarly, as in other reported cases, the system presented here exhibits a 
drastic change in the degree of coherence as a function of the initial particle 
density but in a more abrupt fashion. However, such result does not require the 
use of periodic boundary conditions as it is customary the case in other classes 
of models. 

Further investigation of this model is required in order to develop insights 
on key questions regarding the phenomenon of phase synchronization observed, 
as well as, the dependence of the characteristic time of synchronization on the 
model parameters and on how robustness of synchronization is affected by the 
density of particles. 

On the other hand, since our system can be considered of the minimal ones 
that includes a deterministic steering, it can be readily augmented in order to 
account for more realistic situations. Depending on the context of application 
one might easily consider further the inclusion of particle accelerations, diffcr- 
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cnt types of long or short range interactions for the interparticle attraction, 
inclusion of inhomogeneities and ambient gradients effectuated at the level of 
the interactions and finally the introduction of population variability such as 
groups of leading and/or inert particles. 
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